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CHAPTER  1 


INTRODUCTION 

The  finite  element  aethod  is  a  numerical  analysis  technique  which 
can  be  used  to  find  the  location  in  N- space  at  which  some  scalar 
functional  of  N  variables  is  stationary.  This  aethod  is  one  of  the 
most  popular  methods  used  to  find  solutions  of  both  linear  and  non¬ 
linear  problems  in  the  theory  of  elasticity.  Finite  element  techniques 
which  allow  for  a  systematic  check  on  convergence  with  respect  to  mesh 
size,  discretization  accuracy,  are  useful  to  the  numerical  analyst. 

The  accumulation  of  rcund-cf £  errors  and  the  errors  made  when  numerical 
integration  is  used  can  destroy  the  accuracy  of  a  finite  element 
algorithm.  In  the  case  of  nonlinear  elasticity  the  rcund-off  and 
integration  errors  can  be  magnified  because  of  the  incremental  tech¬ 
niques  used.  When  the  deformations  in  an  elastic  body  are  small  and 
when  the  material  properties  are  not  dependent  on  the  deformation  then 

e 

the  use  of  the  finite  element  method  leads  to  an  energy  functional,  the 
discrete  energy  functional,  which  is  quadratic  in  N-space.  The  first 
variation  of  this  functional  leads  to  a  set  of  linear  equations.  The 
solution  of  these  equations  gives  the  point  in  N-space  where  the 
functional  is  stationary.  If  the  material  properties  are  dependent  on 
the  deformation  or  if  the  deformations  are  large  then  the  discrete 
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energy  functional  is  not  quadratic  in  N-space.  The  first  variation  of 
the  energy  functional  then  leads  to  a  nonlinear  set  of  equations.  The 
solutions  of  these  equations  represent  the  stationary  points  in  N-space 
where  the  discrete  functional  is  stationary.  Since  there  may  be  more 
than  one  solution  the  stability  of  the  discrete  energy  functional  at 
these  points  becomes  important.  Solutions  of  interest,  stable  points, 
are  those  points  in  N-space  for  which  the  discrete  energy  functional  is 
a  local  minimum.  Stability  in  N-space  requires  that  the  Hessian  of  the 
discrete  energy  functional  be  positive  definite  at  a  solution  point. 

1.1  The  Energy  Functional  for  Finite  He  formations  of  Mooney 

Materials. 

A  historical  account  of  the  development  of  the  first  energy 
expressions  used  in  the  analysis  of  rubber  materials  is  given  in 
Treloar's^  book  on  rubber  elasticity.  Due  to  the  fact  that  vulcanized 
rubber  consists  of  long  molecules  linked  together  to  form  an  irregular 
three-dimensional  network  early  researchers  developed  methods  for 
describing  the  behavior  of  a  large  assembly  of  elastic  links.  As 
described  by  Treloar  these  theories  were  developed  during  the  period 
from  1936  to  1947.  The  problem  of  determining  physical  variables  to 
describe  the  energy  functional  was  resolved  when  a  theory  was  developed 
to  describe  the  pure  homogeneous  strain  deformation  of  rubber  by 
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characterizing  the  strain  in  terns  of  three  principal  extension  ratios 
A,  »  ^  *  an<*  ^3  along  three  mutually  perpendicular  axes.  A  brief 

summary  of  Treloar's  description  of  this  theory  is  given  here. 
Definition  (D.l) :  Pure  homogeneous  strain  and  principal  stretch 

ratios  \t  , 

Pure  homogeneous  strain  exists  when  a  unit  cube  in 
the  undeforaed  body  is  transformed  into  a  rect¬ 
angular  parallelepiped  having  three  unequal  edge 
lengths  A.  ,  "\xf  and  1.  See  Figure  1.1.  The 
deformed  edge  lengths  *\(  ,  and  Aj  are 

called  Che  principal  stretch  ratios. 

Definition  0,2") :  Strain  invariants  X  ,  T  ,  and  T  . 

i  ’  —i ’  3 

x,- 

xJ=  V;x\x* 


Assumption  (A.l) :  The  material  is  assumed  to  be  isotropic  in  the 

unstressed  state  and  incompressible  in  all  states 
of  deformation. 

The  requirement  that  the  pure  homogeneous  deformation  be  incompressible 


-  I 


implies  that 


or  that 


-4- 


« 


FIGURE  U  A  PURE  HOMOGENEOUS 

STRAIN  DEFORMATION 
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1=  =  XXV,  - 1 


(1.1.2) 


Then,  for  an  incompressible  material  there  are  only  two  stretch  ratios 


that  are  a  function  of  the  deformation,  X  and  Using  (1.1.2) 


in  (1.1.1)  ve  have 


■J-  v  *- 


X,  '  'K,  +  K  4 


(1. 1.3) 


>:  x 


The  strain  energy  density  functional,  jj  ,  for  a  material  satisfying 


D. 1,  D.2,  and  A.l  is  given  by 

•O  «o 


p  =  Z  I  CVX,'V  (rt-3X 


Cl. 1.4) 


A -o 


where 


C  *  s 


C- 

■i 


are  material  constants. 

The  most  common  form  of  (1.1.4)  used  in  the  analysis  of  finite  deforma' 


tlons  of  incompressible  materials  is  the  Mooney^’ form  which  is 


given  as  follows 
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(J  =  C,[(I 


(1.1.5) 


There  a  re  ocher  formulations  for  finite  deformations  of  incom¬ 
pressible  materials.  Recently  Cescotto  and  Fonder^  have  summarized 
these  formulations.  Their  work  is  concerned  with  variational  prin¬ 
ciples  which  allow  for  a  more  general  behavior  of  the  material  than 
the  Mooney  formulation  allows  for.  Oden  and  Key^^  have  considered 
the  effects  of  using  different  energy  functions  in  finite  elasticity. 
They  noticed  that  by  using  different  strain  energy  functions,  related 
to  the  Mooney  form,  solutions  are  obtained  which  are  significantly 
different.  We  are  concerned  with  the  effect  of  numerical  integration 
and  mesh  refinement  on  the  solutions  to  problems  with  a  given  strain 
energy  formulation.  The  formulation  chosen  here  is  one  of  n-tn-frairr 
potential  energy  using  the  Mooney  form  for  the  strain  energy.  This 
formulation  has  been  shown  to  be  meaningful  by  Tielking  and  Feng^\ 
They  analyzed  nonlinear  .axisyszaetric  membrane  problems  by  applying  the 
Ritz  method  to  the  -total  potential  energy  functional,  TV  ,  given  as 
follows 


(1.1.6) 


where  \j  -■  -the- volume  of  the  undeforaed  elastic  body 

and  VvJ  —  c*le  wor^  ^ne  by  the  applied  fractions. 

We  will  be  concerned  here  with  the  discretisation  of  this  functional 


(4)  (8) 

can  be  found  in  books  by  Oden  ,  and  Strang  and  Fix  .  We  present 
here  a  summary  of  the  finite  element  method  as  it  applies  to  this 
study.  The  technique  presented  for  numerically  evaluating  the  element 
matrices  is  suitable  for  evaluating  the  effects  of  numerical  integra¬ 
tion  on  the  accuracy  of  the  finite  element  solutuons. 

Consider  an  integral  functional  on  (0,  1)  of  the  following  form. 


I(-)  - 


(1.2.1) 


F  may  contain  rational  expressions  of  U.  .uJ  or  their  powers.  Par- 
tition  of  the  interval  (0,  1)  into  N  sub  intervals  (  4{  4 1'+, ) 

where  the  points  }  t  -  l)  2}  •••  }  M+f  are  called  the  nodes. 

Define  a  set  of  variables  •[  u,-  (  which  are  the  finite  element 

approximations  to  U.  and  its  derivatives  at  the  set  of  nodes  ]$,  ]  • 

A* 

That  is,  U.’  is  the  approximation  to  U.  at  node  i  located  at  4;  • 
We  now  wish  to  approximate  U.  between  the  nodes  by  polynomial  inter¬ 
polation  which  is  completely  defined  by  the  nodal  variables  ^  ^  J 


( 


The  power  o£  the  finite  element  nethod  lies  in  the  next  step*  A  low 


order  interpolation,  usually  between  1st  and  5th,  is  chosen  for  the 
interpolation.  This  requires  the  polynomials  to  be  locally  defined 
over  intervals  associated  with  a  small  number  of  nodal  variables  which 
define  the  coefficients  of  the  polynomials*  These  intervals  or  groups 
of  intervals  (  (?)  ,  €  -  f }  •  -  *  )  )  are  called  elements*  When 

this  procedure  is  followed  an  approximation  to  the  integral  (1.2.1)  is 
available  in  the  following  form. 


1(5) 


(1.2.2) 


or 


j  F'Cv)u,u)<)o' 

(?) 


(1*2.3) 


ICS)  is  the  approximation  to  .  The  desired  functions 

U.(/^)  are  those  which  sake  stationary.  The  displacement 

finite  element  solutuons  for  a  given  partition  and  interpolation  order 
are  found  by  using  the  calculus  of  variations  to  determine,  of  all 
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possible  nodal  values  |  Cl^ ■  U ;  j  which  sec  makes  station¬ 

ary*  Thus,  the  variational  boundary  conditions  oust  be  enforced  on 
the  finite  element  approximation.  This  is  usually  a  trivial  task  since 
Che  finite  element  approximation  functions,  as  described  above,  are 
locally  defined  and  modifications  need  be  made  only  over  a  few 
elements*  The  fora  of  the  finite  element  approximation  for  elements 
not  on  the  boundary  remains  unchanged. 

If  the  partition  is  allowed  to  become  successively  finer  then  the 
accuracy  of  the  finite  element  approximations  increases  (until  com¬ 
putational  errors  arise).  The  accuracy  and  rate  of  convergence  of  the 
finite  element  approximations  when  the  partition  is  refined  is  a 
function  of  the  order  of  interpolation  chosen  in  elliptic  boundary 
value  problems*  The  use  of  higher  order  interpolations  (quadratic, 
cubic,  quint ic)  in  linear  elliptic  problem s  often  proves  useful.  In 
the  problems  studied  in  this  dissertation  does  not  have  an 

associated  linear  elliptic  differential  operator  and  the  rates  of  con¬ 
vergence  are  not  known  in  terms  of  parameters  characterizing  a  dif¬ 
ferential  operator* 

Consider  the  discrete  integral  functional  given  by  (1.2.3).  Con¬ 
struct  e  uniform  partition  of  (0,  1)  such  that  ths  intervals  of  ths 


station- 


possible  nodal  values  |  U  4-  j  which  set  makes  J,(u) 
ary.  Thus,  the  variational  boundary  conditions  oust  be  enforced  on 
the  finite  element  approximation*  This  is  usually  a  trivial  task  since 
the  finite  element  approximation  functions,  as  described  above,  are 
locally  defined  and  modifications  need  be  made  only  over  a  few 
elements.  The  fora  of  the  finite  element  approximation  for  elements 
not  on  the  boundary  remains  unchanged. 

If  the  partition  is  allowed  to  become  successively  finer  then  the 
accuracy  of  the  finite  element  approximations  increases  (until  com¬ 
putational  errors  arise).  The  accuracy  and  rate  of  convergence  of  the 
finite  element  approximations  when  the  partition  is  refined  is  a 
function  of  the  order  of  interpolation  chosen  in  elliptic  boundary 
value  problems.  The  use  of  higher  order  interpolations  (quadratic, 
cubic,  quintic)  in  linear  elliptic  problems  often  proves  useful.  In 
the  problems  studied  in  this  dissertation  Xdu.}  does  not  have  an 
associated  linear  elliptic  differential  operator  and  the  rates  of  con¬ 
vergence  are  not  known  in  terms  of  parameters  characterizing  a  dif¬ 
ferential  operator. 

Consider  the  discrete  integral  functional  given  by  (1.2.3).  Con¬ 
struct  s  uniform  partition  of  (0,  1)  such  that  the  intervals  of  the 
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where 


the  variation  of  the  element  nodal  variables. 


z!  {MT{«}em{M’W so 

e*\  L  e 


(1.2.7) 


where 


co 

tai=U{?}t  = 


the  element  gradient  vector 


the  global  gradient  vector 


“■*  w  *  u  m, 


the  global  variation  of  nodal 


variables.  Since  £<£u}  is  arbitrary  except  for  the  constraints 
we  have  the  following  condition  on  { at  a  stationary  locatio 

of  • 


hi  -  hi 


(1.2.8) 


When  finding  the  solution  to  (1.2.8)  and  when  checking  the  stability 
of  the  stationary  point  the  Hessian,  £  J?  3  »  of  is 

.  h...  . 
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often  used*  It  is  given  as  follows* 

r  L' 


(1.2.9) 


The  Hessian  is  computed  on  an  element  bases  as  shown  below 


W  -  U  W 

c  *  i 


U 

0  - , 


*WI 


(1.2.10) 


1.3  The  Use  of  Numerical  Integration  to  Obtain  Element  Matrices. 

Numerical  integration  is  used  in  the  finite  element  method  when 

exact  integration  is  not  practical.  This  is  the  case  when  high  order 

interpolation  is  assumed,  when  multi-dimensional  curved  elements  are 

used,  and  when  the  integral  functional  is  not  a  polynomial  in  the 

(8) 

finite  element  nodal  variables.  Strang  and  Fix  state  the  funda¬ 
mental  problem  as  follows:  "What  degree  of  accuracy  in  the  integra¬ 
tion  formula  is  required  for  convergence?  It  is  not  required  that 

every  polynomial  which  appears  be  integrated  exactly."  This  question 

(9) 

was  answered  for  elliptic  problems  of  2mth  order  by  Fried  in  1974. 
The  requirement  for  elliptic  problems  is  that  the  terms  in  the 
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integral  functional  should  be  integrated  by  a  rule  which  is  exact  for 
polynomials  of  degree  2(p-m)  where  m  is  the  highest  order  derivative 
appearing  in  the  integral  functional  and  p  is  the  order  of  the  inter¬ 
polation  polynomials*  If  we  tried  to  extend  this,  idea  to  the  case 
when  the  Integral  functional  contains  rational  expressions  of  U-fW) 
and  U  (j)  then  we  would  require  exact  integration  to  some  as  yet 
undetermined  order  of  the  rational  expressions*  Such  a  requirement 
would  be  strict  and  nay  lead  to  an  expensive  integration  routine.  Also, 
such  a  requirement  may  not  even  make  sense.  For  example,  the  errors 
involved  in  assuming  the  solution  is  a  polynomial  (the  interpolation 
assumption)  may  be  larger  than  the  errors  introduced  by  inexactly 
integrating  the  rational  expressions  in  the  integral  functional 

A* 

I  L  Iu"\  ]  a  1^  order  scheme.  In  this  study  problems  with 

*  • 

rational  expressions  in  the  energy  integral  are  solved  by  the  finite 
element  method  with  low  order  Integration  schemes.  The  effects  of 
increasing  the  order  of  accuracy  of  the  integration  schemes  are  pre¬ 
sented. 

Let  be  a  mapping  of  the  interval  ^c')  onto 

the  interval  (-1,  1).  Then  the  element  gradient  and  Hessian  matrices 
can  be  numerically  computed  as  follows. 


14 


P 


<i*v 


*FCf,W> 


(1.3.1) 


and 


M,=  I 


iFua-u 


(1.3.2) 


i 


where  j  ^Zj  .  ..  ,  p 

p  •=.  the  number  of  integration  points 
y.  -  the  jth  integration  point  in  (-1,  I) 

■J 

-  the  weight  at  the  jth  integration  point  . 

This  integration  technique  has  recently  been  applied  by  Fried^^  to 
the  nonlinear  finite  element  analysis  of  cantilever  beans  and  plates. 
Element  matrices  containing  rational  expressions  of  the  nodal  variables 
were  evaluated  with  a  low  order  integration  scheme.  The  rank  of  the 
element  Hessian  was  reduced  but  when  the  boundary  conditions  were 
applied  a  full  rank  global  Hessian  was  obtained.  In  this  study  we  will 
determine  numerically  when  the  rank  of  the  element  and  global  Hessians 
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are  reduced. 

1.4  Newton's  Method  and  the  Incremental  Loading 

Fortailation. 

The  method  used  to  determine  the  solution  to  the  nonlinear 
equations  given  by  (1.2.8)  is  a  sequential  method  and  is  called 
Newton's  method.  Given  an  initial  vector  { ulj  the  method  computes 

a  series  of  vectors  {u|  ,  ju ^  which  converge  to  a  vector 

•{  for  which  •  It  is  necessary  that  the 

initial  guess  vector  {  u}^  be  close  enough  to  |u."j  for  a  two 

term  Taylor  expansion  of  to  be  approximately  valid.  For  this 

reason  a  modification  of  the  method,  called  the  incremental  loading 
formulation,  is  used  when  it  is  difficult  to  estimate  a  good  starting 
point  {aI0  •  Newmon's  method  can  be  described  by  considering  a 
Taylor  expansion  of  at  the  Nth  vector  of  the  series.  This  is 

given  as  follows  , 


(1.4.1) 


Requiring  to  be  zero  and  assuming  Is  not 

singular  we  have  Newton's  method  (1.4.2)  for  generating  the  series 
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which  converges  to  { u"J  . 

{«]  ~M  {?] 

1  >Wi  1  V  L  L  N 

When  the  solution  is  dependent  on  a  parameter,  p  ,  which  may  be  large 
it  is  often  possible  to  determine  useful  initial  vectors  {u"^0 
only  when  p  is  small.  In  finite  elasticity  p  is  the  loading. 

If  the  loading  is  large  Newton's  method  is  modified  by  finding  the 
solution  for  a  small  load  P(  and  using  this  solution  as  an  initial 
guess  for  loading  =  P,  t  &P  where  AP  is  a  small  loading 

increment.  This  process  is  continued  until  the  solution  at 
P  -  =  +■  HP  is  obtained.  This  modification  of 

Newton's  method  is  called  an  incremental  loading  formulation. 

Other  methods  can  be  used  to  obtain  solutions  to  (1.2.8).  3oth 
random  searches  and  gradient  techniques  are  popular  alternatives. 

Also,  Newton's  method  as  described  above  can  be  improved^^  but  the 
algorithms  described  here  worked  well  on  the  problems  solved  in  this 

study.  A  more  detailed  discussion  of  techniques  for  solving  nonlinear 

(4) 

equations  in  N-dimensional  spaces  is  available  in  Oden's  book. 


18 


CHAPTER  2 

FINITE  ELEMENT  ANALYSTS 

2.1  The  Potential  Energy  Functional  for  Axisynmetric 

Incompressible  Membranes. 

The  stretch  ratios  for  ax i symmetric  incompressible  membranes 
take  a  form  which  allows  for  a  general  derivation  of  the  finite 
element  matrices  (see  Ref.  7).  Coefficients  in  the  derived  form  of 
the  element  matrices  are  specified  by  the  particular  axisymaetric 
problem  chosen  for  analysis.  Using  the  Mooney  form  (1.1.5)  of  the 
strain  energy  per  unit  volume  with  C(  -  \  and  using  (1.1.6)  we 

have  the  following  fora  of  the  potential  energy  functional. 

TT  (2.i.D 

v 


The  geometry  of  axisyrsaetric  membranes  is  shown  in  Figure  2.1.  We 
assume  the  undefcrmed  membrane  is  of  unit  thickness.  There  is  no 
dependence  on  the  variable  €  shown  in  Figure  2.1  so  (2.1.1)  becomes 


+  oc(irs)]fz Js  -  y 


(2.1.2) 


AXIS  Or 

Symmetry 


( R£)  =  COORDINATES  OF  UNDEFO RMED  SURFAC 
(X,Y).=  COORDINATES  OF  DEFORMED  SURFACE 

FIGURE  2.1  GEOMETRY  OF  AXISYM- 

METRIC  MEMBRANES 


The  stretch  ratios  are  given  as  follows. 

>  AS1  /  (*')" 

A'  *  AS  =  l  i  Cs'^  ) 

(2.1.3) 

-  >.  -  AC1  - 
At  AC  R 


where  /  0\  -  cJt0) 

k  ;  d  R 


The  domain  of  integration  is  now  partitioned  as  described  in  section 
1.2  and  over  each  element  of  the  domain  we  approximate  )(  and  vf 

in  terms  of  the  element  nodal  variables.  Using  the  appropriate  inter¬ 
polation  functions  we  have 


z-mH 


(2.1.4) 


and 


9  =  M». 


Using  these  approximations  the  squares  of  the  stretch  ratios  become 


(2.1.5) 


]h 

A'  "  A 

z  ^ 

where  -  )  +  (S*  )L 


Mapping  the  element  onto  the  interval  (-1,  1)  with  the  relation 
R  =  P(f)  y  using  the  expression  (1.1.3)  for  the  strain 
invariants  X(  and  Xz  ,  and  using  numerical  integration  at 
p  points  as  described  in  section  1.3  we  arrive  at  the  fol 

lowing  form  of  the  discrete  potential  energy. 

Nt 

IT  -  Z  X  (2.1.6) 

e«i 


where 

TTt 


22 


+- 


A,  F;  F; 


<*mri  (uie  )({uC  [  *.r  i  we ) 


HT[!^w,-nwc 


•nd  Wc 


the  term  associated  with  the  work  done  by  the 


applied  tractions 
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2.2  The  Element  Gradient. 

The  expression  for  the  total  potential  energy  can  be  changed  into 
a  form  which  is  more  convenient  for  computing  the  gradient.  The 
stretch  ratios  ^  and  \  ^  in  an  element  are  expressed  in  terms 
of  three  functions  of  ^  .  These  expressions  are  given  below. 

A(M<) 

&(j ,M«)  *  He{vt|>] '  -  B(f ,H)  <2’2"n 

c  (f X  W ■  { wJH a  C  ( i  ,h.) 

and  ~k  ^  become 

(2.2.2) 


Using  (2.2.1)  the  stretch  ratios 

_  /  A1  +  B2  N*7’- 


A 


a,-  t 


The.  total  potential  energy  expressed  in  terms  of  A,  3,  and  C  is 
obtained  by  substituting  (2.2.1)  into  (2.1.6).  vJe  obtain  the  follow- 


T. 


lag  expression  for 


TT  -  Iir  Z  u,,F.F,  (a'  £»»?)  -  F.  c.; 


•  =1 


-3) 


(2.2.3) 


*^a7f:vy*‘  <-s*k‘  <■  a,f-7/)Ys‘  )' 


4-  F.'c"1  -  3 

»  I 


-  w, 


The  element  gradient  is  now  found  by  differentiating  (2.2.3)  with 
respect  to  {u\e  •  have 

.  (v,„l 

)  _  r  tpctjT 

^w;  "*  ,J 


Using  (2.2.4)  we  obtain,  after  simplification,  the  element  gradient. 
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We  =  jh; 

^  i  it 


/ 


(2.2.5) 


{AiHl*  e.WKl^PVKA;  F^M  F •>>»*,)  c7 

.  C(fo]  C I  -  *A,  A!F.-F>^»l)'crjS 


5W« 

iWTt 


2.3  The  Element  Hessian 

The  expression  developed  for  the  element  gradient  can  be  used  to 
compute  the  element  Hessian.  We  again  use  the  expressions  (2.2.1)  to 

obtain 


c  ^Me) 

Mc 


*8(T,M«) 

$ 


(2.3.1) 
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and 


aCq.M,) 

3He 


■  {ffp5 


Substi&ittion  of  (2.3.1)  into  (2.2.5),  differentiation,  and  sisal if Ica- 
tion  leads  to  the  elenent  Hessian. 


Ms 

t*I.  - 


(2.3.2) 


0**^0 (.A.'  C-  ^  p  ■  0>  <Tc?(l  (A.)  *  s  • 


*(l+«.F.C))(A,-F  -A.F,  F.  c;Y'  (A’1**)  ')[{'V.!i4'.\ 


tT 

j 


4  V <  A; f: VaV ») ))(F,\jA.. f)  F,’/A>  6*  )'C)  [  ] J  <ft] 


+<(C;k;  f;v.1 


*Z8.'(W»s‘)cVA..F.tF;.V<<C.A: 
+‘FA;6ic7(A>6))A;FtF.l(i  +  « 


jWe 
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2.4  Quadratic  and  Cubic  Elements. 

The  particular  forms  of  interpolation  within  an  element  can  affect 
the  results  of  an  analysis.  The  quadratic  three-node  element  and  the 
cubic  two-node  elements  were  used  in  this  study.  The  quadratic  three- 
node  element  has  only  displacements  for  nodal  variables  and  is  a  C° 
element.  The  cubic  two-node  element  has  both  displacements  and  first 
derivatives  as  nodal  variables  and  is  a  element.  These  one-dimen¬ 
sional  elements  have  interpolation  functions  related  to  Lagrange  and 
Hermite  polynomials.  The  quadratic  element  is  called  a  second  order 
Lagrange  element  and  the  cubic  element  is  called  a  third  order  Semite 
element. 

The  quadratic  three-node  element  has  the  following  set  of  nodal 
variables,  see  Figure  2.2. 


(2.4.1) 


The  associated  interpolation  functions  on  (-1,  1)  are 

{fvnf 1  ff/n, 0  >V's>  >  °  ^ 

(2.4.2) 
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f= element  coordinate 

R= GLOBAL  COORDINATE 


FIGURE  2.2  QUADRATIC  AND  CUBIC 

ELEMENTS 
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T 

and  (4ty]  °  .^H),  ° 

'f/A'i  -  -j_  ^  n  -  o 

4^0--  i-T 

•Pj/f)  ■  iK'-f) 


The  cubic  two-node  element  nodal  variables  are 


The  associated  interpolation  functions  on  (-1,  1)  are 


and 

W<vfr{°>° 


e>,  o,%H) 


%f«r 


\X.  - 


(2.4.3) 


(2.4.4) 


fhere 


•P/o  "  i(i  -  MM  ) 


*ftfo 


V" 


--  M'+nV-O 


em  *  fft-OC'-f.^ 


and  /\  -  I/j_ 


the  length  of  the  element  in  the  global 


system. 


2.5  •  Boundary  Conditions. 

When  stationary  values  of  the  integral  functional  given  rn 
(1.2.1)  are  to  be  determined  it  is  necessary  to  enforce  the  essential 
boundary  conditions  vhich  are  determined  using  variational  calculus. 
The  essential  boundary  conditions  for  (1.2.1)  are 


=■ 


(2.5.1) 


at  the  boundaries.  If  U.  is  not  specified  at  a  boundary  then 

*  O  there.  In  the  finite  element  analysis  this  require¬ 
ment  determines  a  set  of  equations  vhich  the  nodal  variables  oust 


satisfy.  The  requirement  that  the  gradient, 


» 


be  zero  still 


applies.  A  gradient  equation  with  a  reduced  dimension  rust  then  be 
found.  The  new  equation  specifies  the  location  of  stationary  values 
of  the  discrete  energy  functional  for  which  the  nodal  variables  also 
satisfy  the  boundary  conditions. 

The  axisymnetric  problems  analyzed  in  this  study  are  one  dimen¬ 
sional  boundary  value  problems.  The  boundary  conditions  for  such 
problems  specify  relations  among  the  nodal  variables  in  only  the  first 
and  last  elements.  Then,  the  computation  of  the  gradient  and  Hessian 
is  modified  only  for  these  boundary  elements.  Each  boundary  equation 
is  arranged  to  express  the  nodal  variable  at  the  boundary  in  terms  of 
the  remaining,  reduced,  element  nodal  variables.  When  this  has  been 
done  we  can  determine  the  full  set  of  nodal  variables  in  terms  of  the 
reduced  nodal  variables.  We  have 

HtE[c]eH,«  <2-5-2> 


where 


the  reduced  set  of  nodal  variables 


and 


With  expression  (2.5.2) 


the  matrix  of  coefficients  determined  by  the 
constraint . 

computation  of  the  reduced  global  gradient  and 
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Hessian  is  easily  accomplished.  If  we  substitute  (2.5.2)  into  (2.2.1) 


we  arrive  at 

-K. I ■M’W.W.; 


(2.5.3) 


C  ( i  ,Ht«)  W„  r  cUrt>}'{  ?m}  [  c  ]JuJ -■  C  ( f ,  i  uy 

CX  t  C  ^ 


The  element  gradient  and  Hessian  will  then  have  the  same  scalar  coef¬ 
ficients  as  those  given  in  (2.2.5)  and  (2.2.2).  Only  the  vectors  and 
matrices  need  to  be  modified  and  we  arrive  at 


h.  ■  i c  i;  w« 

W„'(C]'[*],[CJ< 


(2.5.**) 


CHAPTER  3 


AX  I  SYMMETRIC  TEFCRMATICK  0?  A  DIEK 
WITH  A  HOLE 

3.1  Introduction. 

(12) 

In  1951  Rivlin  and  Thomas  studied  the  in-plane  behavior  of  a 
thin  disk  of  Mooney  material.  The  disk  contained  a  central  circular 
hole  and  was  radially  loaded  on  its  outer  edge.  Rivlin  and  Thomas* 
analysis  compared  well  with  experimental  test  data.  They  determined 
that  a  value  of  0.1  for  0(  in  equation  (1.1.5)  leads  to  computa¬ 
tional  results  that  are  in  good  agreement  with  experimental  test  data. 
Using  the  cylindrical  coordinate  systems  shown  in  Figure  3.1  Rivlin 
and  Thomas  assumed  that  the  sheet  was  thin  enough  so  that  the  deforma¬ 
tions  are  not  dependent  on  2.  Further,  they  considered  only  axisym- 
metric  deformations,  which  eliminates  dependence  on  0  .  The  result 

is  that  the  axisymmetric  deformation  of  a  thin  disk  is  a  pure  strain 
deformation  with  the  stretch  ratios  dependent  on  only  R.  They 
obtained  an  approximate  analytical  solution  to  this  finite  deformation 
problem  as  follows.  Expressions  for  all  the  radial  and  tangential 
stretch  ratios  ‘A,  and  \z  were  determined  in  terms  of  the  inde¬ 
pendent  coordinate  ,  the  dependent  coordinate  R,  and  the  form  of 

the  internal  energy  W.  These  expressions  were  used  with  the 
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D  =  ENFORCED  DISPLACEMENT  OF  OUTER  EDGE 
Hr  THICKNESS  OF  SHEET  =  1.0 
R  =  UNDEFORMED  . COORDINATES 
X  =  DEFORMED  COORDINATES 


FIGURE  3.1  DISK  WITH  HOLE. 
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equilibrium  equation  and  boundary  conditions  to  determine  the  first 
three  coefficients  of  a  Taylor  expansion  of  the  solution.  The  Taylor 


expansion  being 


f>(R) 


4. 


N  •»  o 


(3.1.1) 


where  ^  =  the  internal  radius  of  the  disk. 

A  numerical  solution  to  this  problem  was  found  by  assuming  knowledge  of 
and  “X  at  a  particular  radius  R  =  R^-.  The  equilibrium 
equation  and  a  differential  equation  obtained  by  eliminating  the  inde¬ 
pendent  variable  from  the  expressions  for  "X  t  and  "X  ^ 

were  used  with  the  incompressibility  condition  XT"X  ~  I 

to  obtain  approximate  expressions  for  Xt  and  ~XZ  at  nearby 

values  of  R  and  thus  for  all  values  of  R. 

(13) 

Recently  Veraa  and  Rana  considered  the  axisymmetric  in-plane 
deformation  of  a  disk  with  a  hole.  They  used  a  strain  energy  function 
which  is  different  than  the  Mooney  form.  The  energy  function  used  is 


w  -  t-ie  -  BX, 


(3.1.2) 


J 
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where 


I*  = 


X  +■  \  ^  A-s  -  3 


N 


and  G>  S  —  material  constants. 

When  N  =  |  and  fA  -  1  (3.1.2)  reduces  to  (1.1.5)  with 

c(  -  0.6  .  Under  these  conditions  on  N  and  M  Verma 

and  Sana  reduce  their  general  equation  to  a  special  equation  called  the 
Varga  equation.  They  compare  the  solution  of  this  equation  to  both 
Rivlin  and  Thomas'  three-term  Taylor  expansion  of  the  solution  for 
Mooney  material  and  experimental  data.  The  three-term  expansion  com¬ 
pared  better  with  experimental  data.  This  may  be  because  using 

N  -  t  and  ]  in  (3.1.2)  implies  a  slightly  different 

form  of  the  strain  energy  than  the  Mooney  form.  In  the  Mooney  form 

the  value  of  C(  was  chosen  to  make  the  computed  and  measured  data 

« 

agree  well. 

Finite  element  solutions  to  the  axisyrsnetric  deformation  of  an 

incompressible  disk  with  a  hole  have  been  obtained  by  a  different 

finite  element  technique  than  presented  in  this  study.  In  1967 
(14) 

Oden  presented  a  method  for  obtaining  numerical  solutions  to  prob¬ 
lems  in  the  theory  of  nonlinear  elasticity.  He  used  the  finite  element 
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method  to  formulate  the  nonlinear  problem  and  solve  it.  Most  of  the 
finite  element  formulations  used  previous  to  that  of  Oden  were  based 
on  finding  successive  corrections  tc  the  linear  problem.  Oden  discret¬ 
ized  the  domain  of  the  elastic  body  into  subdomains  called  finite 
elements.  The  discretization  was  assumed  fine  enough  so  that  the  dis¬ 
placements  vary  linearly  within  an  element.  The  expressions  describing 
the  linear  interpolation  of  the  displacements  inside  an  element  in 
terms  of  the  nodal  displacements  were  used  to  obtain  expressions  for 
the  strains  in  an  element.  These  strain  expressions  are  used  with  the 
nonlinear  theory  of  elasticity  to  obtain  the  strain  energy  in  an 
element  in  terms  of  the  nodal  displacements.  Expressions  for  the 
generalized  forces  acting  cn  an  element  are  derived  and  the  total 
potential  energy  for  the  finite  element  is  computed.  The  condition  of 
minimum  potential  energy  is  used  to  determine  a  set  of  nonlinear  equi¬ 
librium  equations  for  the  element.  A  global  coordinate  system  is  used 
and  the  appropriate  transformations  are  applied  so  that  the  element 
equilibrium  equations  are  expressed  in  the  global  system.  The 
assembled  global  nonlinear  equations  are  then  determined  and  solved. 

A  solution  of  the  disk  problem  analyzed  in  this  chapter  by  Oden's 
method  is  given  in  reference  4,  chapter  IV. 
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3.2  Analysis  of  a  Disk  with  a  Hole. 

In  this  section  we  analyte  the  disk  problem  stated  in  section  3.1. 
The  disk  has  an  inner  radius  of  1.0  and  an  outer  radius  of  3.0.  The 
boundary  conditions  are  an  enforced  displacement  cf  the  outer  edge  and 
a  stress-free  inner  edge.  The  finite  element  method  outlined  in  chap¬ 
ter  2  is  used.  The  thickness  of  the  disk  and  the  material  constant 
in  (1.1.5)  will  both  be  unity.  The  outer  edge  of  the  disk  will  be 
given  a  prescribed  displacement  and  an  initial  guess  of  the  solution 
will  be  used  to  start  the  Newton  algorithm.  The  quadratic  three-noce 
element  will  be  used  with  two  and  four  point  Gauss- Legendre  integra¬ 
tion. 


Referring  to  Figure  3.1 
the  tangential  stretch  ratio. 


the  radial 


stretch  ratio, 
can  be  computed. 


Tv ,  »  and 

They  are 


~  dR 


(3.2.1) 


Using  the  incompressibility  condition  (1.1.2)  the  third  stretch  ratio 

is 


(3.2.2) 

i 

i 

i 

The  radius,  R,  is  the  only  independent  variable  and  X  is  the  only 
dependent  variable.  Partition  the  domain  of  R  into  a  uniform  mesh  for 
quadratic  three-node  elements  each  of  length  2.  (  .  Each 

element  is  mapped  to  ^  coordinates  belonging  to  the  interval 
(-1,  1)  by  the  mapping  function 


* 


R  =  uc  4- 


(3.2.3) 


where  Y{c  =  the  radial  displacement  to  the  center  node  of 

the  undefonned  element. 

These  relations  were  used  in  the  finite  element  analysis  described  in 
chapter  2  to  compute  the  numerical  data  described  in  this  chapter.  To 
start  the  Newton  algorithm  all  nodes  were  moved  radi  ally  outward  by 
the  amount  the  outer  edge  was  being  displaced.  With  no  nodal  force 
associated  with  the  innermost  node  the  traction-free  boundary  con¬ 


dition  at  the  inner  radius  is  satisfied.  The  outermost  element's  nodal 


displacements  must  be  reduced  in  number  Because  of  the  enforced  em¬ 
placement  at  the  outer  edge.  If  the  enforced  displacement  of  the 
outermost  node  is  defined  in  the  initial  guess  then  the  outermost 
element's  gradient  and  Hessian  can  be  modified  as  follows  and  the 
enforced  displacement  boundary  condition  will  be  satisfied. 


{?}  =  [  c  ]  [9 ' 

L  Nfc  L  JN  [ 


(3.2.4) 


where 


[c] 


m  the  outermost  element's  identification  number 


M. 


)  O  o 


O  I  o 


The  finite  element  program  used  to  obtain  the  data  for  this  prob 
lem  was  written  in  double  precision.  The  Newton  algorithm  was 
terminated  when  the  ^2.  norm  of  the  gradient  was  less  than  10~^. 
The  norm  of  the  displacement  increments  was  less  than  10~^ 
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when  the  algorithm  was  terminated  indicating  that  the  displacement 
vector  had  converged  in  the  Newton  algorithm  to  more  than  eight  sig¬ 
nificant  digits. 

The  data  presented  by  Rivlin  and  Thomas  in  reference  12  was  used 
to  define  four  problems.  The.  stretch  ratios,  7>-.  ,  at  R  =  3.0 

were  read  from  the  four  curves  on  the  graph  of  vs  R  in  refer¬ 
ence  12.  These  values  of  were  used  to  determine  four  values 

for  the  outer  edge  disp 1 acement,  D,  of  the  disk,  shown  is  Figure  3.1. 

The  curves  of  vs  R  which  result  from  the  finite  element 

analyses  are  shown  in  Figure  3.2.  The  changes  of  7\z  vs  R  with 
mesh  size  and  integration  order  are  not  detectable  graphically  and 
these  differences  are  discussed  in  the  next  section.  The  agreement 
between  Rivlin  and  Thomas'  computations  and  the  finite  element  com¬ 
putations  for  E00^*  Tiie  distribution  of  the  strain 

invariants  T  and  I  as  predicted  by  the  finite  element 

« 

analyses  is  shown  in  Figures  3.3  and  3.4.  These  curves  also  agree 
well  with  the  data  presented  in  reference  12. 

3.3  Convergence  of  the  Finite  Element  Solutions. 

To  study  the  convergence  rate  of  the  displacement  approximations 
the  rate  of  convergence  of  the  inner  radius  of  the  disk  with  respect  to 


*-v  . 
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a  uniform  nesh  refinement  was  determined.  Ve  as  suae  that  we  have  the 
following  power  fora  of  convergence. 


✓V* 


(3.3.1) 


where 


UN  ~ 
U  a 


P  - 


C  = 


the  finite  element  inner  radius  for  N  elements 

the  converged  finite  element  inner  radius 

the  element  sire  for  N  elements 

the  convergence  rate 

a  constant  which  is  not  dependent  on  the  mesh 
size. 


Vhen  (3.3.1)  is  used  with  the  data  in  Table  3.1  the  convergence  rate, 
p  ,  and  the  converged  displacement,  u_  ,  can  be  determined  for 
each  integration  order.  The  converged  inner  radius  for  two-point 
integration  was  equal  to  the  converged  inner  radius  for  four-point 
integration  to  seven  digits  (see  Table  3.1).  The  rank  of  the  element 
Hessians  was  determined  in  the  first  and  last  iterations.  They  were 
full  rank  in  both  cases.  Also,  the  converged  solution  had  a  positive 
definite  Hessian.  Figure  3.5  shows  the  convergence  rates  and  accuracy 


TABLE  3.1  Convergence  Data  for  Disk  with  Hole 


Inner  Radius  (D  = 

5.10) 

No.  of  Integration 

Points 

6  Elements 

9  Elements 

12  Elements 

Converged 

2 

5.9543663 

5.9543443 

5.9543404 

5.9543385 

4 

5.9543421 

5.9543393 

5. 954338S 

5.9543386 

Strain  Enerev  (D  = 

5.10) 

No.  of  Integration 
Points 

6  Elements  9  Elements  12  Elements  Converged 

2 


302.08010  302.08307 


302.08359 


302.08384 


LOG 


log10(ne) 


FIGURE  3.5  CONVERGENCE  OF  INNER 

RADIUS-DISK  WITH  HOLE 


of  the  finite,  element  solutions  associated  with  using  both  two-point 
and  four-point  integration  for  computing  the  element  matrices.  The 
convergence  rates  are  nearly  identical  but  four-point  integration  gave 
better  accuracy  by  a  full  order  of  magnitude. 

The  convergence  rate  of  the  energy  approximations  was  determined 
by  the  sane  method  that  was  used  to  compute  the  convergence  rate  for 
the  displacement  approximations.  Figure  3.6  shows  that  the  convergence 
rate  using  two-point  integration  was  only  slightly  greater  than  the 
convergence  rate  associated  with  four-point  integration.  However,  the 
accuracy  of  the  energy  approximations  associated  with  four-point 
integration  was  an  order  of  magnitude  better  than  the  accuracy 
associated  with  two-point  integration. 


LOG 


LOG0(N£) 

0.0  0.5  1.0  1.5 


FIGURE  3.6  CONVERGENCE  OF  STRAIN 

ENERGY-DISK  WITH  HOLE 


CHAPTER  4 

INFLATION  cf  A  circular  disk 
4.1  Introduction. 

The  inflation  of  an  initially  flat  circular  disk,  see  Figure  4.1, 

C 15 ) 

made  of  incompressible  material  was  analyzed  by  Adkins  and  Rivlin  A 
and  Green  and  Adkins  It  is  assumed  that  the  deformed  disk  has 

principal  radii  of  curvature  that  are  everywhere  large  compared  with 
the  thickness  of  the  disk.  This  allows  the  stress  variation  over  the 
thickness  of  the  disk  to  be  neglected  in  the  analysis.  Adkins  and 
Rivlin  used  the  six  differencial  equations  derived  by  Love^7'  to 
represent  the  equilibrium  of  a  thin  sheet.  Four  of  these  equations  are 
identic lly  satisfied  when  the  symmetry  conditions  are  introduced 
leaving  only  two  equations  to  be  solved.  The  two  principal  tractions 
are  the  dependent  variables,  and  the  deformed  radius  is  the  independent 
variable.  Expressions  for  these ^tractions  are  derived  in  terms  of  the 
stretch  ratios  and  the  strain  energy  function.  The  stretch  ratios  are 
then  determined  in  terms  of  the  undeformed  and  deformed  coordinates. 

An  approximate  solution  is  then  found  by  a  technique  similar  to  the  ' 
technique  used  to  find  the  approximate  solutions  for  the  disk  with  a 


hole  mentioned  in  Chapter  3 
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A  different  analysis  of  the  axisymmecric  deformations  of  nonlinear 

(13) 

membranes  was  made  by  Yang  and  reng  *  in  1970.  In  their  paper  the 
stretch  ratios  "X  and  X  ^  in  the  meridian  and  circumferential 
directions  respectively  are  parametrically  expressed  in  terms  of 
coordinates  and  ^  (R)  ,  and  a  ^function  ^  which 

specifies  the  undeformed  configuration.  The  parametric  variable  ^ 
is  the  radial  coordinate  from  the  axis  of  symmetry.  The  equilibrium 
equations  of  an  axisymmetric  membrane  loaded  with  an  internal  pressure 
and  a  constitutive  equation  expressing  the  relationship  between  the 
principal  tractions  end  stretch  ratios  are  used.  Three  new  functional 
variables  U.  ,  AT  ,  and  uu~  are  defined  in  terms  of  2  ,  r;  , 


7% 


,  and 


/ 0 


so  that  the  ecuations  of  equilibrium  reduce  to 


a  set  of  first  order  ordinary  differential  equations.  In  the  analysis 
of  the  inflation  of  a  circular  disk  they  are  able  to  reduce  the  problem 
to  a  nondimens ional  form  which  allows  the  solutions  to  be  obtained  if  a 
stretch  ratio  at  the  pole  is  assumed.  That  is,  if  a  stretch  ratio  at 
the  pole  is  assumed  then  a  configuration  and  pressure  is  determined 
which  produces  the  assumed  stretch  ratio  at  the  pole.  Their  results  • 
compared  well  with  the  earlier  work  of  Adkins  and  Rivlin^^X  Tielking 

and  Feng*'')  used  the  potential  energy  principle  and  the  Kits  method  in 
1974  to  determine  the  deformation  profiles  of  an  inflated  disk  and  thei: 
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results  compare  well  with  Yang  and  Feng's  work. 

The  finite  inflation  of  a  circular  disk  has  been  analyzed  by  the 
finite  element  method.  In  reference  4  Oden  gives  a  summary  of  the  work 
done  from  1966  to  1971.  In  Oden's  analysis  the  displacement  interpola¬ 
tion  was  linear  and  the  strain  energy  function  involved  the  strain 
invariants  in  an  exponential  form.  The  exponential  form  of  the  strain 
energy  used  approximates  the  Mooney  form  for  the  material  constants 
used  in  the  analysis.  The  results  of  the  finite  element  analysis  were 
compared  with  experimental  data  and  the  agreement  was  good.  Oden  noted 
the  highly  nonlinear  character  of  this  problem  by  presenting  a  graph  cf 
the  inflation  pressure  vs  the  polar  extension  ratio.  Certain  values  of 
pressure  have  more  than  one  associated  polar  extension  ratio  indicating 

that  the  deformation  is  not  uniquely  defined  in  terms  of  the  polar 

(19) 

extension  ratio.  Recently  Argyris  et  al  considered  the  use  of 
higher  order  elements  for  the  analysis  of  the  finite  inflation  of  a 
circular  disk  made  of  Mooney  material.  The  data  obtained  using  higher 
order  elements,  quadratic  interpolation,  was  compared  to  data  obtained 
using  linear  interpolation.  The  results  indicated  that  the  use  of  quad¬ 
ratic  interpolation  elements  leads  to  better  accuracy  with  fewer 
elements. 
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4.2  Analysis  of  the  Inflation  of  a  Circular  Disk. 

The  inflation  of  an  initially  flat  circular  disk  of  unit  radius 
is  analyzed  here.  The  Mooney  fort  of  the  energy  is  assumed  with 
C. .  -  l  .  C>  and  ^  i  o.  1  in  equation  (1.1.5).  Due  to  the 
symmetry  involved  only  the  radial  domain  (o,  1)  has  to  be  discretized. 

The  boundary  conditions  are  zero  displacement  at  the  outer  edge  and  zero 
slope  at  the  pole.  The  incremental  pressure  -  formulation  is  used.  The. 
solution  algorithm  was  started  with  an  initial  out-of-plane  deformation. 
The  quadratic  three-node  element  is  used  with  two-  and  four-point  Gsuss- 
Legendre  integration. 

We  use  X.  as  the  radial  displacement  and  Y  as  vertical 
displacement  cf  the  points  along  the  radius  of  the  undeformed  disk,  see 
Figure  4.1.  Then,  the  meridional  stretch  ratio,  A  ,  ,  and  the  circum¬ 
ferential  stretch  ratio,  7i  ^  »  can  computed.  They  are 


(4.2.1) 


I 


A 


2_ 


and 


R 
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The  incompressibility  condition  t*A  ^  = 

stretch  ratio 


yields  the  third 


R 

x( 

The  radius,  R  ,  is  the  only  independent  variable.  ^  and  Y 
are  the  functions  of  R  to  be  approximated.  The  domain  of  R.  is 
partf.tior.ed  for  a  uniform  mesh  of  quadratic  three-node  elements  as  in 
chapter  3.  The  mapping  function  (3.2.3)  is  used  and  the  finite  element 
solutions  are  determined  by  the  method  described  in  chapter  2.  The  con¬ 
straint  equations  which  represent  the  boundary  condition  at  the  pole, 

R  =  q  ,  can  be  computed  with  the  aid  of  Figure  4.2.  In  element 
coordinates  we  have 


vtv  --  jwnl'W, 


Y(p*4^r-oy,  +o-f >viy(VO  <4-2-*r 
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we  require 


(4.2.5) 


and 


Using  (4.2.4)  and  (4.2.5)  we  obtain 


(4.2.6) 


or 


(4.2.7) 


The  constraint  relation  for 


the  element  at  the  outer  edge  of  the  disk  is 
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The  discrete  fora  of  (4.2.10)  in  element  coordinates  with  p  inte¬ 
gration  points  is  given  as  follows. 

P 

We  -  C;  8;  (4.2.11) 

» •  l 

where  C  and  "0  are  given  in  (2.2.1). 

The  element  gradient  and  Hessian  work  expressions  calculated  from 
(4.2.11)  are 


(4.2.12) 


and 


o>  We 


~  -2irf 


The  stretch  ratios  ft  and  A  ^  evaluated  at  the  pole  are. 
equal  since  the  circumferential  stretch  ratio  becomes  a  meridional 
stretch  ratio  at  chat  point.  The  stretch  ratio  "/(  at  the  pole  can 
not  be  computed  directly  from  (4.2.1)  and  the  nodal  displacement  data 
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since  the  radius  ^  is  zero  there.  In  element  coordinates  we  have 


for  the  element  at  the  pole*  The  stretch  ratio,  "X  ^  ,  at  the  pole 

then  becomes 


n  1  V 

rs  Ai- 


-  -i-  X3 

26  3 


(4.2.13) 


Equations  (4.2.1)  to  (4.2*13)  were  used  in  the  finite  element 

analysis  described  in  chapter  2  to  compute  deformations  of  the  inflated 

disk.  The  computer  program  was  written  in  double  precision.  Cross- 

secticnal  profiles  of  the  disk  at  three  pressures  are  shown  in  Figure 

« 

4.3.  These  profiles  were  drawn  using  displacement  data  associated  with 
9  elements.  The  resulting  stretch  ratios  for  the  same  pressures  are 
shown  in  Figures  4.4  and  4.5.  The  stretch  ratios  were  drawn  from  data 
computed  using  3  elements.  Nodal  data  was  graphed  and  a  smooth  curve 
fit  to  the  nodal  data  over  each  element.  The  stretch  ratios  associated 
with  9  elements  at  J  -  S,  S*0  are  also  shown  to  graphically  present 


FIGURE  4.4 


A,  VS  R  FOR  IN 
FLATED  DISK. 
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the  changes  computed  in  the  stretch  ratios  as  a  result  of  mesh  refine¬ 
ment.  The  figures  indicate  that  a  mesh  near  9  elements  yields  an 
accuracy  sufficient  for  most  applications.  Also,  it  appears  that  for  a 
coarse  mesh  the  accuracy  is  best  near  the  first  and  last  nodes  of  these 

three-node  one-dimensional  elements.  The  distributions  of  the  strain 

«• 

invariants  for  £  «  S. S’O  are  shown  in  Figure  4.6. 

4*3  Convergence  of  the  Finite  Element  Solutions. 

The  rate  of  convergence  of  the  displacement  and  energy  approxima¬ 
tions  with  respect  to  a  uniform  mesh  refinement  was  determined  by  the 
same  method  used  in  section  3.3.  The  computed  pole  displacement  and 
total  strain  energy  values  at  a  pressure  p  -  S.SO  and  at 
several  mesh  values  are  shown  in  Table  4.1.  The  computer  program 
written  to  obtain  the  data  in  Table  4.1  was  written  in  double  precision, 
The  computed  converged  data  for  two-  and  four-point  integration  in  Table 
4.1  compares  well  to  six  significant  digits  for  the  pole  displacement 
and  to  five  significant  digits  for  the  strain  energy.  The  pole  dis¬ 
placement  is  not  a  finite  element  approximation  variable  but  is  a 
linear  combination  of  two  approximation  variables  through  equation  • 
(4.2.6).  The  results  of  the  convergence  calculations  are  displayed  in 
Figures  4.7  and  4.8.  The  higher  order  integration  was  not  helpful  in 


FIGURE  4.6  STRAIN  INVARIANTS  VS 
.  R  FOR  INFLATED  DISK. 


TABLE  4.1  Convergence  Data  for  Inflation  of  Disk 


No.  of  Integration  t 

Points  3  Elements  6  Elements  9  Elements  Converged 


1.7012760  1.7321020  1.735031 5  1.7355636 

1.6862435  1.7321128  1,7348283  1.7355563 


26.278508  26.402439  26.408690  26.410114 


901 
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FIGURE  4.7  CONVERGENCE  OF  POLE 


DISPLACEMENT  -  INFLATED 
DISK. 


LOG 


FIGURE  4.8  CONVERGENCE  OF  STRAIN 

ENERGY -INFLATED  DISK. 
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improving  hhe  accuracy  cf  die  computed  crspiacement  rield  and  in  race, 
gave  slightly  less  accurate  results.  The  sane  results  apply  to  the 
accuracy  cf  the  strain  energy.  However,  the  C--;- -point  integration 
scheme  ;as  an  order  of  magnitude  more  accurate  than  the  ic-ur-point 

integration  scheme  for  computation  of  the  tonal  strain  energy. 

U.b  Out -of -Plane  Deformation  cf  a  Unit  Disk  with  a  Concentric 

Rigid  Inclusion. 

The  in-plane  axi symmetric  deformations  of  a  rubber  sheet  with  a 

(1S») 

circular  Tiigbd  incYvisbcn  bp.s  bc^rn  invost/igHtsd  by  Y  in  Yf'b T 

and  again  by  Verna  and  Rana^-^)  in  1??S.  The  out-cf-plane  deforma¬ 
tions  due  to  a  vertical  axial  force  on  the  rigid  inclusion  and  a 
pressure  on  the  membrane  were  determined  ir.  137h  by  Tic  Using  and 
Feng^b  In  this  section  the  cut -of -plane  deformations  resulting 
from  a  prescribed  pole  displacement  of  the  rigid  inclusion  are  deter¬ 
mined,  see  Figure  b.c.  3oth  a  large  and  a  small  inclusion  were 
considered . 

The  computer  program  used  to  obtain  the  data  for  the  inflation  o 
a  unit  disk  was  modified  to  allow  the  domain  of  R  to  be  C  0 
&3  shown  in  Figure  h.9.  The  initial  guess  of  the  solution  specified 
the  (X,  T)  location  of  the  cuter  edge  of  the  inclusion.  The  first 
element  had  its  gradient  and  Hessian  contributions  reduced  with  the 
following  relations  to  enforce  the  boundary  condition  at  the  outer 
edge  of  the  rigid  inclusion. 
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FIGURE  4.10  CROSS  SECTION  PROFILES 

FOR  DISK  W\TW  WCLOSIOU. 
INCLUSION  RADIOS  =  0.3. 
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FIGURE  4.12.  STRETCH  RATIOS  FOR  DISK 

WITH  INCLUSION,  INCLUSION 
RADIUS  *  0.3  . 
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CHAPTER  i> 

INFLATION  OF  A  TCRUS 
5.1  Introduction. 

It  is  difficult  to  obtain  analytical  solutions  valid  for  the 
finite  inflation  of  a  torus  (see  Figure  5.1)  because  a  torus  has  two 
principal  curvatures.  The  earliest  analyses  of  the  inflated  torus 
involved  solving  stress  equilibrium  equations  derived  in  the  linear 

theory  of  elasticity.  These  solutions  result  in  a  discontinuous  dis- 

(20'' 

placement  field  when  the  stress  field  is  integrated.  In  1963  Jordan  ' 
reformulated  the  torus  problem  by  allowing  snail  deformations  but  large 
derivatives  of  the  deformations.  Instead  of  approaching  only  the 
stress  equations  Jordan  added  the  requirement  that  the  displacement 
field  be  continuous.  The  material  constitutive  law  was  Hooke’s  law  for 
an  isotropic  material  in  two  dimensions.  The  formulation  included 

i 

defining  a  special  "primary  shape  function"  which  allowed  him  to  write 
the  equilibrium  equation  and  continuity  equation  in  a  form  which  could 
be  solved  by  an  iterative  procedure.  This  shape  function  was  used  in 
later  years  by  other  investigators  to  aid  them  in  their  computations  of 
approximate  solutions  for  the  nonlinear  inflation  (large  strains)  of  a 
torus.  Jordan's  solution  was  the  first  solution  to  have  a  continuous 
displacement  field.  He  found  that  his  solutions  for  the  hoop  stress 
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FIGURE  5.1  TORUS  WITH  CIRCULAR 

CROSS  SECTION. 


72 


agree  well  with  the  earlier  solutions  but  that  the  overall  torus  cir¬ 
cumferential  stress  predicted  by  earlier  solutions  was  in  error  by  up 
to  187.. 

(21) 

In  1963  Sanders  and  Liepins  used  a  nonlinear  membrane  theory 
(small  strains  but  large  rotations)  to  formulate  the  nonlinear  equi¬ 
librium  equations  for  a  toroidal  membrane  under  internal  pressure.  The 
material  constitutive  law  used  is  Hooke's  law  for  an  isotropic  material 
in  two  dimensions.  Sanders  and  Liepins  were  able  to  convert  the  non¬ 
linear  equilibrium  equations  inco  a  form  which  could  be  solved  by 
asymptotic  methods.  Their  solutions  agree  well  with  Jordan' s^*^ 

solutions  and  their  technique  was  easier  to  apply  than  Jordan's. 

(22) 

In  1965  Liepins  determined  the  natural  vibration  modes  of 
prestressed  toroidal  membranes.  The  stresses  resulting  from  the 
linear  analysis  of  toroidal  membranes  were  used  to  define  the  prestress 
state.  The  equations  for  the  analysis  of  vibrations  of  prestressed 
membranes  were  reduced  to  a  set  of  second  order  differential  equations 
by  separation  of  variables.  The  separated  equations  represent  different 
forms  of  vibration  of  the  torus;  that  is,  flexural  modes  due  to  the 

torus  bending,  extensional  modes  of  the  meridional  curve,  etc. 

(23) 

Also  in  1965  Kydoniefs  and  Spencer  analyzed  the  finite  infla¬ 
tion  of  an  incompressible  elastic  torus.  They  stated  the  theory  for  the 
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uniform  inflation  of  a  thick  wall  torus.  Solutions  were  obtained  for 
the  case  when  the  cross  sectional  radius  of  the  torus  is  very  small 
when  compared  with  the  overall  torus  radius.  They  presented  seme 
numerical  results  for  the  case  of  neo-Eookean  material.  Kydoniefs  and 
Spencer  suggest  that  since  the  number  of  exact  solutions  to  problems 
in  finite  elasticity  is  small  "and  seems  likely  to  remain  so"  that 
approximate  solution  techniques  for  these  problems  are  of  interest. 
Their  interest  was  with  the  use  of  perturbation  procedures.  In  par¬ 
ticular  they  studied  series  expansions  of  the  solution  in  terms  of  a 
numerically  small  geometric  parameter  (the  ratio  of  the  torus  cross 
section  radius  to  the  overall  radius  of  the  torus).  Their  solutions 
were  valid  only  for  the  case  when  the  cross  section  remains  circular 

and  is  small  compared  to  the  overall  radius  of  the  torus. 

(24) 

In  1967  Kydoniefs  and  Spencer  computed  asymptotic  solutions 

for  the  finite  inflation  of  a  toroidal  membrane  of  an  initially 

circular  cross  section  and  made  of  a  Mooney  material;  see  equation 

(1.1.5).  The  parameter  in  the  asymptotic  expansion  was  the  ratio  of 

the  radius  of  the  circular  cross  section  to  the  overall  radius  of  the 

torus.  As  a  result,  their  solutions  apply  only  to  the  case  when  the 

torus  has  a  small  cross  sectional  radius.  They  used  the  "shape 

(20) 

function"  introduced  by  Jordan  to  help  simplify  the  algebra. 


The  zero,  first,  and  second  order  equations  resulting  from  the 
asymptotic  expansions  are  solved  in  their  paper.  In  describing  the 
solutions  for  a  torus  cade  of  Mooney  material  Kydoniefs  and  Spencer 
indicate  that  there  exists  a  pressure  at  which  the  volume  of  the  torus 
will  increase  without  bound  with  no  further  increase  in  pressure.  In 
this  chapter  the  finite  element  solutions  indicate  the  same  type  of 
response  for  a  torus  with  a  large  cross  sectional  radius  compared  to 
its  overall  radius. 

5.2  Analysis  of  the  Inflation  of  a  Torus  of  Circular  Cross 

Section. 

The  deformations  resulting  from  the  inflation  of  a  torus  which 
has  a  circular  cross  section  in  its  undeformed  state  (see  Figure  5.1) 
are  analyzed  in  this  section.  The  Mooney  form  of  the  strain  energy 
as  given  in  equation  (1.1.5)  is  used  with  C{  -  |.C  and  o(  *  0.1 
The  cross  section  radius  a  a  0.5  and  the  overall  radius  of  the  torus 
R  »  1.0.  The  rectangular  caordinates  of  the  unde formed  torus  are 

(x,  y)  and  of  the  deformed  torus  are  (X,  Y)  as  shown  in  Figure  5.1. 
Using  the  angular  coordinate  co  as  the  independent  variable  and  con¬ 
sidering  the  symmetry  we  determine  that  the  domain  O  £  wo  ± 
must  be  discretized.  The  boundary  conditions  are 
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Y  (to  =  o 


and 


cU 

ci  uo 


-  O 

Ui  =  OsTT 


(5.2.1) 


Using  the  parametric  coordinates  q  and  ^  to  measure  the  arc 
length  of  the  undeforaed  and  deformed  meridional  curve  respectively, 
the  meridional  stretch  ratio,  "h  (  ,  and  the  circumferential  stretch 

ratio,  ,  can  be  computed.  They  are  given  by  the  following 

expressions 


>  ■  JL 

1  V 
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the  expressions  for  A,  and  in 


terms  of  X,  Y  and  uj  become 

(X")1  4-  (Y')1 


^  - 


(R  +  a,  Co-au})*" 


(5.2.3) 
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and 


When  we  use  the  quadratic  elements  the  boundary  conditions  expressed 
by  equations  (5.2.1)  yield  the  following  relations  among  the  nodal 
variables  for  the  first  and  last  elements  on  the  domain  O  i  uj  6  Tr 
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The  first  and  last  element  contributions  to  the  gradient  and  Hessian 
are  appropriately  modified  with  relations  (5.2.4)  to  (5.2.7)  as  shown 
in  Chapter  2. 


The  work,  Uj  ,  done  by  the  pressure,  p  ,  during  the  de forma¬ 


tion  is 


U 


-  j  f<)V  =  4-ttP^  X  Y  -  ?V. 


(5.2.8) 


U>s*Tr 


where  V  =  the  initial  volume  of  the  torus.  Crooning  the  const; 
tern  since  it  does  not  contribute  to  the  gradient  or  Hessian  and 
charging  the  integration  limits  we  have 


M 

r 

UJ  -  H-irf  \  XT*  J 


(5.2.3) 


ui  =  o 


Using  J  ( 4  X  )  -  XX  q  Ui  ,  integrating  by  parts,  and  remembering 
that  'Y(Ui=0i‘n')  =  O 


we  obtain 


TT 
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(5.2.10) 
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(5.2.11) 


where  W&  *  the  contribution  to  the  work  from  element  (e). 

Expression  (5.2.11)  is  identical  in  form  to  (4.2.10)  and  after  some 
algebra  ve  obtain  the  following  element  contributions  to  the  gradient 
and  the  Hessian 


(5.2.12) 


and 
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1  *e  1  «,,  L  V  J 


(5.2.13) 


Ihe  expressions  for  the  stretch  ratios  X,  ,  ,  and  /l»3 

In  equations  (5.2.3)  were  used  with  the  boundary  conditions  shown  in 
(5*2.4)  to  (5.2. 7),  and  the  contributions  from  the  work  shown  ix 
(5.2.12)  and  (5.2.13)  in  the  finite  element  algorithm  outlined  in 
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Chapter  2  to  compute  the  deformation  of  the  torus  as  a  function  of 
pressure.  The  computer  program  vas  written  in  double  precision.  The 
cross  sectional  profiles  of  the  inflated  torus  are  shown  in  Figure  5.2. 
It  was  found  that  the  inner  circumferential  radius  of  the  torus 
decreases  slowly  with  increasing  pressure  while  the  outer  circumferen¬ 
tial  radius  increases  rapidly  with  increasing  pressure,  see  Figure  5.3. 
The  distribution  of  the  stretch  ratios  and  strain  invariants  for  three 

values  of  the  pressure  are  shown  in  Figures  5.4  and  5.5.  In  the  paper 

(24) 

by  Kydoniefs  and  Spencer  the  following  expressions  for  the  membrane 
stresses  are  presented  (here  the  initial  membrane  thickness,  hiQ  , 
is  assumed  to  be  1.0). 


t,  ■  xn 


(5.2.14) 


ax. 


The  stress  distributions  for  three  values  of  the  pressure  were  com¬ 
puted  using  (5.2.14)  and  are  shown  in  Figure  5.6. 

The  maximum  values  of  T(  and  were  computed  for  a 

series  of  values  of  pressure  and  the  results  are  shown  in  Figure  5.7. 


ai 


AXIS  OF  TORUS 


FIGURE  5.2  INFLATION  OF  TORUS. 


PRESSURE  P 


FIGURE  5.3  INNER  AND  OUTER  TORUS 

RADII  VS  PRESSURE. 
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FIGURE  5.5  STRAIN  INVARIANTS  FOR 
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FIGURE  5.6  STRESS  RESULTANTS 

FOR  INFLATION  OF 
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FIGURE  5.7  INFLATION  OF  TORUS, 

MAX  MEMBRANE  STRESS 
VS  PRESSURE. 
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When  the  pressure  increased  beyond  £  s  3.50  it  became  difficult 

to  use  the  incremental  pressure  formulation  to  find  configurations  at 

successively  higher  pressures.  As  a  result,  the  computations  were 

stopped  at  £  =  3.5  23  .  This  behavior  was  suggested  by 

(24) 

Kydoniefs  and  Spencer  when  they  determined  asymptotic  solutions 
for  the  inflation  of  a  torus  with  a  small  cross  section  and  made  of 
Mooney  material. 

5.3  Convergence  of  the  Finite  Element  Solutions. 

The  convergence  of  the  displacements,  strain  energy,  and  the 
lowest  eigenvalue  of  the  Hessian  of  the  potential  energy  are  examined 
in  this  section.  As  in  section  3.3  convergence  was  assumed  to  follow  a 
power  law  with  respect  to  mesh  size  as  shown  in  equation  (3.3.1).  The 
outer  overall  circumferential  radius  of  the  torus  and  the  total  strain 
energy  were  determined  for  a  series  of  mesh  sizes.  The  quadratic 
element  was  used  with  2 -  and  4-point  integration  and  the  cubic  element 
w as  used  with  3-,  4-,  end  6-point  integration  (2-point  integration 
with  the  cubic  elements  resulted  in  a  reduced  rank  Hessian).  The 
finite  element  programs  written  to  obtain  the  data  used  in  this  section 
were  written  in  double  precision.  Table  5.1  shows  the  convergence  of 
the  outer  radius  and  the  strain  energy  when  quadratic  elements  are  used. 
The  results  indicated  that  both  2-  and  4-point  Gauas-Legendre 


TABLE  5.1  Convergence  Data  for  Inflation  of  Torus,  Quadratic 


Elements 

Outer  Radius  (?  a  3.25) 


No.  of  Integration 
Points 

.3  Elements 

6  Elements 

9  Elements 

Converged 

2 

1.7538039 

1.7591435 

1.2 593895 

1.7594454 

4 

1.7541697 

1.7591665 

1.7593966 

1.7554489 

i 

i 


Strain  Er.ergv  (? 

a  3.25) 

No.  of  Integration 
Points 

3  Elements 

6  Elements 

9  Elements 

Converged 

2 

12.589861 

12.689202 

12.692323 

12.692851 

4 

12.687924 

12.695023 

12.693598 
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integration  gave  comparable  accuracy  vhen  the  quadratic  element  was 
used.  Table  5.2  indicates  similar  results  vhen  cubic  elements  vere 
used.'  In  fact,  to  eight  significant  digits  there  were  no  differences 
in  the  computed  displacements  for  the  cases  of  4-  and  6-point 
integration  vhen  the  cubic  elements  vere  used.  The  converged  values 
of  the  outer  radius  and  strain  energy  vere  identical  to  eight  signifi¬ 
cant  digits  for  3-,  4-,  and  6-point  integration  with  cubic  elements. 
Thus,  the  data  in  Tables  5.1  and  5.2  indicates  that  using  an  integra¬ 
tion  scheme  of  higher  order  than  required  for  a  full  rank  Hessian  is 
wasteful  on  these  problems. 

Figures  5.8  and  5.9  summarize  the  accuracy  of  the  integration 
schemes  used.  These  figures  are  not  intended  for  comparing  the 
efficiency  of  the  cubic  element  vs  the  quadratic  elements  since  the 
total  degrees  of  freedom  are  not  accounted  for  in  these  figures. 
However,  comparable  convergence  rates  vere  obtained  for  both  the 
quadratic  and  cubic  elements.'  This  is  surprising  since  in  linear 
elliptic  problems  this  convergence  rate  is  dependent  on  the  order  of 
the  interpolation  polynomials  (see  reference  8). 

The  eigenvalues  of  the  Hessian  vere  investigated  and  the  data  is 
presented  in  Table  5.3.  The  Hessian  vas  positive  definite  at  the 
solution  f  U  «j)|  <  \  S ^  |  for  all  mesh  sizes  and  for  all 
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TABLE  5.2  Convergence  Data  for  Inflation  of  Torus,  Cubic  Elements 


Outer  Radius  (?  =  3.25) 


No.  of  Integration 

Points  4  Elements 

8  Elements 

12  Elements 

Converged 

3 

1.7596120 

1.7594601 

1.7594520 

1.7594500 

4 

1.7596131 

1.7594602 

1.7594520 

1.7594500 

6 

1.7596131 

1.7594602 

1.7594520 

1.7594500 

Strain  Enerzv  (? 

*  3.25) 

1 

No.  of  Integration 

Points 

4  Elements 

8  Elements 

12  Elements 

Converged 

3 

12.692719 

12.693058 

12.693070 

12.693072 

4 

12.691811 

12.693046 

12.693069 

12.693072 

6 

12. '691827 

12.693046 

12.~693069 

12.693072 
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TABLE  5.3 


« 

Convergence  Data  for  Inflation  of  Tom,  Lowest  Eigen- 
va'ue  of  Sessian 


Quadratic  Elements  (P  «  3*25) 


No.  of  Integration 
Points 

3  Elements 

6  Elements 

9  Elements 

Converged 

2 

15.969049 

7.3660452 

4.8047576 

0JD507 

4 

16.060752 

7.3717412 

4.805522 5 

0.6935 

Cubic  Elements  (? 

-  3^25) 

*1 

No.'  of  Integration 

•  Pul:. is 

4  Elements 

8  Elements 

12  Elements 

Converged 

3 

3.3165205 

1.6808296 

1.1262479 

-0.0249 

4 

3.3219768 

1.6825444 

1.1269123 

-0.0267 

6 

3.3209281 

1.6824904 

1.1269044 

•0.0262 

integration  orders.  However,  if  the  eigenvalues  are  assumed  to  be 
converging  according  to  a  power  law  with  respect  to  mesh  size,  only 
the  quadratic  element  is  associated  with  a  positive  definite  converged 
solution. 

5.I4.  Analysis  of  the  Inflation  of  a  Torus  of  Elliptical  Cross 

Section. 

This  section  describes  a  finite  element  analysis  of  the  inflation 
of  a  torus  with  an  elliptical  cross  section  in  its  undeforced,  state. 

The  torus  is  assumed  to  have  its  center  at  X  *  R,  senimajer  axis  in  the 
x  direction  of  value  "a",  and  a  semiminor  axis  in  the  y  direction  of 
value  "b"  (see  Figures  5*1  and  5*10).  Following  the  analysis  of 
section  5*2  end  taking  the  geometry  of  the  ellipse  into  account  in  the 
calculations  the  following  form  of  the  stretch  ratios  are  obtained. 

.1) 


V 


[  *  *  f  ]L 


(5.fc.2) 
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FIGURE  5.10  INFLATION  OF  A  TORUS 

WITH  AN  ELLIPTICAL 
CROSS  SECTION!. 
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where 


I 


a.'  V>** 


The  third  stretch  ratio  is  given  "by  the  condition  that  ^  A  1^.3  *  1 
The  expression  used  for  computing  the  work  done  "by  the  pressure  when 
a  torus  with  a  circular  cross  section  is  inflated  applies  here  also. 

0 

This  is  true  because  the  only  difference  in  phe  work  expressions  is  a 
constant  tern  which  does  not  contribute  to  the  gradient  or  the  Hessian. 

The  finite  element  program  written  to  analyte  the  torus  with  a 
circular  cross  section  using  cubic  elements  was  modified  using  expres¬ 
sions  (5*^.1)  and  (5.U.2).  Another  minor  modification  of  the  program 
.was  required  to  represent  integration  on  a  torus  with  an  elliptical 
cross  section  instead  of  a  circular  cross  section.  The  case  of  a 
semimajcr  axis  equal  to  0.5  and  a  semiminor  axis  equal  to  0.3  was 
analyzed.  The  cross  section  profiles  are  shown  in  Figure  5*10.  The 
elliptical  cross  section  very  quickly  converted  to  a  near  circular 
cross  section  when  the  pressure  was  increased  from  0.0  to  about  1.35* 
In  fact,  at  low  pressures  ^  P  <  |.3SQ  the  solutions  obtained 
could  not  he  considered  to  apply  to  a  thin  rubber  membrane  since  a 
compressive  membrane  stress  was  computed  in  the  circumferential  direc¬ 
tion  on  the  outermost  portions  of  the  torus. 

After  a  pressure  of  1.35  was  reached  many  of  the  characteristics 
of  the  deformation  were  similar  to  the  case  of  a  torus  with  a  circular 
cross  section.  Figure  5*11  Indicates  that  the  innermost  and  outermost 
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FIGURE  5.11  INNER  ANI>  OUTER  TORUS 

RAt.ll  VS  PRESS  ORE  FOR 
ELLIPTICAL  CROSS  SECTION. 


circumferential  radii  tend  to  initially  change  magnitude  in  the 
opposite  sense  from  the  case  for  a  circular  cross  section.  Eventually, 
however,  these  radii  follow  the  same  type  of  pattern  as  shown  in 
Figure  5.3.  When  the  pressure  was  large  enough  >  2. 09  )  to 

cause  >  1. 10  everywhere  the  cross  section  profile  was  nearly 
circular  and  Figures  5 -12  to  5 *15  indicate  that  the  stretch  ratios, 
strain  invariants  and  stresses  have  profiles  that  are  similar  to  the 
corresponding  profiles  for  the  torus  with  a  circular  cross  section 
shown  in  Figures  5*^  to  5-7. 
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